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Abstract 

The problem of efficiently finding the best match for a query in a given set with respect to the 
Euchdean distance or the cosine similarity has been extensively studied in literature. However, a closely 
related problem of efficiently finding the best match with respect to the inner product has never been 
explored in the general setting to the best of our knowledge. In this paper we consider this general 
problem and contrast it with the existing best-match algorithms. First, we propose a general branch- 
and-bound algorithm using a tree data structure. Subsequently, we present a dual-tree algorithm for 
the case where there are multiple queries. Finally we present a new data structure for increasing the 
efficiency of the dual-tree algorithm. These branch-and-bound algorithms involve novel bounds suited 
for the purpose of best-matching with inner products. We evaluate our proposed algorithms on a variety 
of data sets from various applications, and exhibit up to five orders of magnitude improvement in query 
time over the naive search technique. 

1 Introduction 

In this paper, we consider the problem of efficiently finding the best-match for a query from a given set of 
points with respect to the inner-product similarity. Formally, we consider the following problem: 
Problem. For a given set of N points S C and a query q G R®, efficiently find a point p E S such that: 

= max((?,r). (1) 

res 

We call this the problem of maximum inner-product search. The focus of this paper is to improve the 
efficiency of this search. An alternate formulation of the above problem in terms of a vector and matrix 
multiplication is as following: 

Problem. For a given vector w G and a matrix M e R®^^, efficiently compute the following: 

||?«^M||_^ = niax(u;'^M). (2) 

This problem appears to be very similar much existing work in literature. Efficiently finding the best 
match with respect to the Euclidean (or more generally Lp) distance is the widely studied problem of fast 
nearest- neighbor search in metric spaces [S] . Efficient retrieval of the best match with respect to the cosine 
similarity is the extensively researched in the field of text mining and information retrieval jTj. But as 
we will explain in the next section, the maximum inner-product search is not only different from these 
aforementioned tasks, but also arguably harder. 



1.1 Applications 

An obvious application of maximum inner-product search stems out of the widely successful matrix- factorization 
framework in recommender system challenges like the "Netfiix prize" [HI HU [5] . The matrix- factorization 
task obtains accurate representation of the available data in terms of user vectors and items vectors (ex- 
amples for items would be movies or music). In this setting, the preference of a user for an item is the 
inner-product between the corresponding user's vector and the item's vector. The efficient retrieval of rec- 
ommendations for a user is equivalent to the problem in equation [1] with the user as the query and the 
items as the reference set. For the challenges, linear scan of the items are usually employed to find the 
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best recommendations. An efRcient search algorithm would make the retrieval of recommendations in the 
matrix-factorization framework scalable to real world systems. 

The usual document retrieval tasks use the cosine-similarity to match documents. However, in certain 
setting [TT], the documents are represented as (not necessarily normalized) vectors and the inner-product 
between these vectors represent the similarity between the documents. In this case, unless the vectors are 
normalized to have the same length, document matching using the cosine-similarity jTj might make the 
algorithm scalable at the cost of returning inaccurate solutions since the inner-product is not the same as 
the cosine-similarity (we will explain this more elaborately in the section [2]). 

There is a similar problem known as the the max-kernel operation: for a given set of points S and a 
query q and a kernel-function /C(-, •), the task is to find the point p G S with the maximum value of K{q,p) 
over the set S. This problem is widely used in maximum-a-posteriori inference |20| in machine learning, and 
for the task of image matching |23| in computer vision. If the kernel function can be explicitly represented 
in the form a function i^(-) such that JC{q,p) = {(f{q),ip{r)), then this problem reduces to the problem in 
equation [T] after all the points in the set S and the query q is transformed into the i^j-space. 



In this paper, we consider the general problem of efficient maximum inner-product search and propose 
two tree-based branch-and-bound algorithms along with a new data structure to solve this problem more 
efficiently that the naive linear scan. In the following section, we contrast this problem to the usual problems 
of nearest-neighbor search in metric spaces and best-matches with respect to cosine-similarity. This presents 
the need for explicit attention to this problem (equation [T]). However, we do motivate the use of the existing 
tree data structures for solving this task efficiently. In section [31 we propose a simple branch-and-bound 
algorithm using the existing ball-tree data structure and a novel bound. In the following section (section 
S]), we address the situation where there are multiple queries on the same set of points and propose a dual- 
tree branch-and-bound algorithm along with a new tree data structure, cone trees, to index the queries. 
The proposed algorithms are evaluated for their efficiency over a variety of data sets in section [H Section |6] 
demonstrates how the proposed algorithms can be applied to the max-kernel operation with a general kernel 
function where it is not required to have an explicit representation of the points in the (yS-space. In the final 
section (section [7]) , we provide our conclusions along with possible future directions for this work. 

2 Maximum Inner-product Search 

The inner-product between two vectors is very closely related to the Euclidean distance between the points 
represented by these vectors as well as to the cosine-similarity between these to vectors. Numerous techniques 
exists for nearest- neighbor search in Euclidean metric space (see surveys like [9]). Large scale best matching 
algorithms have also been developed for the cosine-similarity measure [1], with a lot of focus on text data. 
The problem of nearest-neighbor search (in metric space) has also been solved approximately with the widely 
popular Locality- sensitive hashing (LSH) method |14U18] . The LSH technique has also been extended to other 
forms of similarity functions (as opposed to the distance as a dissimilarity function) like the cosine similarity 
[7P . The approximate max-kernel operations can also be solved efficiently with LSH under certain conditions 
on the kernel function. Some other techniques like dimension reduction |30| and dual-tree algorithms [20\ 
have also been used to solve the approximate max-kernel operation efficiently. 

2.1 How is maximum inner-product search different from existing problems? 

In what follows, we will try to show that the problem stated in equation [T] is different from these existing 
problems. Hence techniques applied to these problems (like LSH) cannot be directly applied to this problem. 

^An important thing to note here is that the similarity function used in Charikar et.al.[7] is not exactly the cosine-similarity. 
The distance between two points p and q was measured by 6 /it, where 8 is the angle made the two points at the origin, making 

the similarity function ( 1 — — ) . This similarity function has a direct correspondence to the cosine similarity. 



1.2 This Paper 
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Figure 1: Best matches: For a given query q, pc, Pe and pi denote the best match with respect to the 
Cosine-similarity, the Euchdean distance and the Inner-product respectively. It is apparent from this figure 
that even on a plane, the best match with respect to these similarity functions can be very different. 



Nearest-neighbor Search in Euclidean Space. The problem of finding the nearest-neighbor in this 
setting can be posed as finding a point p Cz S for a query q such that: 



p = arg min 1 1 g — r 1 1 = arg max 
re5 - r-es 



argmax(g, r) (unless \\r"'^ 

reS 




Hence, if the norms of all the points in S are normalized to have the same length, then the problem of finding 
the best match with respect to the inner-product is equivalent to the problem of finding the nearest-neighbor 
in Euclidean metric space. However, without this restriction, the two problems can have potentially very 
different answers (figure [2]). 

Best-matching with Cosine-similarity. The problem of finding the best match with respect to the 
cosine-similarity can be posed as finding a point p £ S for a query q such that 

{q, t) {q, r) 

p — arg max t]— 777]— 77 — arg max ■ 



res \\q\\ \\r\\ ^ res \\r\\ 

^ argmax(g,r) (unless ||r|| = ky r Cz S). 
res 

As in the previous case, the best match with cosine similarity is the best match with inner-products if all 
the points in the set S are normalized to have the same length. Under general conditions, the best matches 
with these two similarity functions can be very different (see figure [2]). 

Locality-sensitive Hashing. LSH has been applied to a wide variety of similarity functions. LSH involves 
constructing hashing functions such that each hash function h must satisfy the following for any pair of points 
r,p e S: 

Pr[/i(r) = h{p)] = sim(r,p), (3) 
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where siin(r,p) G [0, 1] is the similarity function of interest. For our situation, we can normalize our data set 
such that V r S 5, ||r|| < lH, and assume that the all the data is in the first quadrant (so that none of the 
inner-products go below zero). In that case, sim(r, p) ~ {r,p) e [0, 1] is a valid similarity function of interest. 

It is known that for any similarity function to admit a locality sensitive hash function family (as defined 
in equation 131) , the distance function d{r,p) 1 — sim(r,p) must satisfy the triangle inequality (Lemma 
1 in [7]). However, the distance function d{r,p) — 1 — (r,p) does not satisfy the triangle inequality (even 
when all the points are restricted to the first quadrant So LSH cannot be applied to the inner product 
similarity function even when we assume that all the data lies in the first quadrant (which is quite a restrictive 
assumption). 

Efficient Max-kernel Operation. There are different existing techniques of solving this problem effi- 
ciently. For kernel functions with very high (possibly infinite) dimensional explicit representations, Rahimi, 
et.al., 2007 [3D], propose a technique to transform these high-dimensional representations into lower- dimensions 
while still approximately preserving the inner-product to improve scalability. However, the final search still 
involves a linear scan over the set of points for the maximum inner-product or a fast nearest-neighbor search 
under the assumption that finding the nearest-neighbor is equivalent to maximizing the inner-product. For 
translation invariant kernelfQ, a tree-based recursive algorithm has been shown to scale to large sets [2IJ| . 
However, it is not clear how this algorithm can be extended to the general class of kernels. LSH is widely 
used for image matching in computer vision [2 3) , but only for kernel functions that admit a locality sensitive 
hashing function [7]. 

Hence, none of the existing techniques can be directly applied to our problem (equation [1]) without 
introducing inaccurate results or limiting assumptions. 

2.2 Why is maximum inner-product search possibly harder? 

Unlike the distance functions in metric space, inner products do not induce any form of triangle inequality 
(even under some assumptions as mentioned in the previous section). Moreover, this lack of any induced 
triangle inequality causes the similarity function induced by the inner products to have no admissible family 
of locality sensitive hashing functions. And any modification to the similarity function to conform to widely 
used similarity functions (like Euclidean distance or Cosine-similarity) will create inaccurate results. 

Moreover, inner-products lack a very basic property of generally used similarity functions - the self 
similarity is high (generally the highest). For example, the Euclidean distance of a point to itself is 0; the 
cosine-similarity of a point to itself is 1. The inner-product of a point x to itself is ||a;|| , which may be high 
or low depending on the value of the ||x||. Moreover, there can possibly be many other points like y in the 
set such that (y,x) > 

Hence, without any assumptions, this problem of obtaining the best match with respect to the maximum 
inner product is inherently harder than the previously dealt similar problems. This is possibly the reason 
why there is no existing work for this problem without any restrictions on the domain (at least to the best 
of our knowledge). 

2.3 Are trees the answer? 

In this paper, we explore the tree data structure for indexing the points and a branch-and-bound algorithm 
specifically for the task of maximum inner-product search. Tree data structures have been widely used for 
the task of nearest- neighbor search [ISjIllIlQ]. And even though the task of nearest- neighbor search is slightly 
different from the task of maximum inner-product search, we believe that trees can still be useful for this 
task. 

^This normalization is different than the normahzation mentioned before where all the points were normalized to have the 
same length. Here the lengths are normalized to be less than equal to one, but not equal to each other. 

^Consider the following counter example: Let x,y,z G 5 be points such that [|a;|| = \\y\\ = \\z\\ = 1, and angles made 
between x k. y, y h z and z h x at the origin are (I -0.1), f and (f -0.3) respectively. In this case the inequality, 
d(a;, y) + d(y, z) > d(z, x) does not hold for ■) = 1 - (■, ■). d(a;, y) = 0.226, d{y, z) = 0.293 & d{z, x) = 0.704. 

^Kernel functions K{j>, q) which are dependent only the (Euclidean) distance between the points p and q are considered 
translation invariant kernels. The Gaussian RBF kernel is such a translation invariant kernel function. 



4 



Trees are known to be good indexing schemes in low to medium dimensions, while some new tree data 
structures have been developed for data in high dimensions with some low dimensional structure jlOi |4] . A 
hierarchical representation of the data is useful. In this paper, we try to solve the problem of exact maximum 
inner-product search. The hierarchical tree data structure provides a very intuitive extension to solve the 
problem approximately to gain efficiency [H . 

Moreover, if the search is to be performed with strict constraints - error constraints or time constraints, 
the tree-based branch-and-bound algorithms can be easily adapted for that purpose. This is because these 
branch-and-bound algorithms are incremental algorithms. This is not possible with something like LSH 
- LSH provides theoretical error bounds, but there is no way of ensuring the error constraint during the 
search. Moreover, LSH is inherently not an incremental algorithm, and hence cannot be used in a limited 
time setting. 

An important advantage of trees is that the trees require a single construction - the branch-and-bound 
algorithm adapts for the different levels of approximate and/or time limitations. Hashing techniques require 
multiple hashes for different levels of approximation. The usual norm is to pre-hash for multiple values of 
approximation. Trees can also be constructed by learning from the data using techniques from machine 
learning I^S] to provide better accuracy and efficiency. 

This is why we use trees to solve the problem. Trees might not be the best possible way to solve this 
problem, but trees do bring a lot of advantages with them. 

3 Tree-based Search 

Ball trees [25] are binary space-partitioning trees that have been widely used for the task of indexing 
data sets. Every node in the tree represents a set of points and each node is subsequently indexed with a 
center and a ball enclosing all the points in the node. The set of point at a node is divided into two disjoint 
sets which form the child nodes. This partitions the space into (possibly overlapping) hyper-spheres. The 
tree is built hierarchically and a node is declared to be a leaf node if it contains a set of points of size below 
a threshold value Nq. 




Figure 2: Ball-trees: All the points are limited within the root ball (the bold-face circle). However, the 
subsequent balls does not necessarily lie within the parent ball - the points still lie within the root ball, but 
the ball enclosing the points in the child node are not necessarily compact enough to lie within the parent 
ball. However, the child node would be confined within the parent node if we used hyper-rectangles instead 
of balls to index the data. 

3.1 Tree Construction 

We use a simple ball tree construction heuristic that approximately picks a pair of pivot points which are 
farthest apart from each other [28 , and splits the data by assigning the points to their closest pivot. The 
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intuition behind this heuristic is that these two points might he in the principal direction. The sphtting and 
the recursive tree construction algorithm is presented in Algorithms [T] & [2] for completeness. 

The tree is very space efficient since every node only stores the indices of the item vectors instead of 
the item vectors themselves. Hence the matrix for the items is never duplicated. Another implementation 
optimization is that the vectors in the items' matrix are sorted in place (during the tree construction) such 
that all the items in the same leaf node are arranged serially in the matrix. This is to avoid any random 
access to the memory when accessing items in the same leaf node. 



Algorithm 1 MakcBaUTrccSplit(Data S) 
Pick a random point x G 5* 
A ^ argmax^/gs ||x - 

B ^ argmaxx'gs 11-4 - x'Hj 
return {A,B) 



Algorithm 2 MakeBaUTree(Set of items S) 
Input - Set S 
Output - Tree T 
T.S ^ S 
T.^ mean (5') 
T.R ^ maxpgs \\p - T.fj.\\l 
if 15"! < No then 

/ / Leaf node 

return T 
else 

// else split the set 

iA,B) ^ MakeBallTreeSplit(S') 

Si ^ {pes: \\p~A\\l<\\p-B\\l} 

Sr^S\Sl 

T.lc ^ MakeBallTree(5'i) 
T.rc MakeBallTree(S'r) 
return T 
end if 



Figure 3: Ball-tree Construction: The object T.S denotes the set of points in the node T. T.fj, denotes 
the Euclidean mean of the items in the node T and T.R denotes the minimum radius of the ball centered 
around T./i enclosing all the points in the node T. T.lc and T.rc denotes the left and right child of the tree 
node T. 



3.2 Branch-and-bound algorithm 

Ball trees are widely used for the task of nearest neighbor search and are known to be fairly scalable 
to moderately high dimensions [281 126) . The search usually employs the depth- first branch-and-bound 
algorithm. A nearest neighbor query is answered by traversing the tree in a depth-first manner by first going 
down the node closer to the query and bounding the minimum possible distance to the other branch with the 
triangle-inequality. If this bound is greater than the current neighbor candidate for the query, the branch is 
removed from computation. 

An analogous greedy depth-first algorithm can be used for maximum inner-product search. But instead 
of traversing down the node closer to the query, the choice is made on the basis of the maximum possible 
inner-product between the query and any potential point from the node. The recursive depth-first branch 
and bound algorithm is presented in Algorithm H) The search algorithm for a query (q) begins at the root 
of the tree (Alg. [5]). At each step, the algorithm is at a tree node (T). It checks if the maximum possible 
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inner-product between the query and any point in the node, MIP(g, T), is any better than the current best- 
match for the query (g.bm). If the check fails, this branch of the tree is not explored any more. Otherwise, 
the algorithm recursively traverses the tree, exploring the branch with the better potential candidates in a 
depth-first manner. If the node is a leaf, the algorithm just finds the best-match within the leaf with a linear 
search (Alg. [3]). This algorithm ensures that the exact solution (i.e., the best-match) is returned by the end 
of the algorithm. 



Algorithm 3 LinearSearch(Query q, Reference Set S) 
for each p £ S do 
if {q,p} > q.\ then 
q.hm p 
q.\ <r- {q,p) 
end if 
end for 



Algorithm 4 TreeSearch(Query g, Tree Node T) 

if q.\ < MIP(g,r) then 

// This node has potential 
if isLeaf{T) then 

LmearSearch(g, T.S) 
else 

// best depth first traversal 

I, ^ MIP(<3,r.lc); Ir ^ MIP(g,r.rc); 

if Ii < Ir then 

TreeSearch(g, T.rc); TreeSearch(g, T.lc); 
else 

TreeSearch(g, T.lc); TreeSearch(g, T.rc); 
end if 
end if 
end if 

// Else the node is pruned from computation 
return; 



Algorithm 5 FindExactMaxIP(Query set V, Reference Set S) 

T ^ MakeBallTree(5) 
for each q £V do 

q.X i c«; 

q.hm <— 0; 

TreeSearch(g, T); 

return q.hm; 
end for 



Figure 4: Single-tree Search: The object q.hm contains the current best-match candidate for the query 
and q.X denotes the inner-product between the query q and its current best-match q.hm. The function 
'MIP(g, T) = {q, T.fi) + \\q\\ T.K denotes the upper bound on the maximum possible inner-product between 
the query q and any point lying in the tree node T. 
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Figure 5; Bounding with a ball 
3.2.1 Bounding maximum inner-product with a ball 

Since the triangle inequality does not hold for the inner product, we present an novel analytical upper bound 
for the maximum possible inner product of a given point (in this case, the query q) with points in a ball. It 
is important to note that the information about the ball is limited to its center and its radius. For the rest 
of this section, we use the notation ||-|| to denote the ||-||2. 

Theorem 3.1. Given a ball Bpg of points centered at pq with radius Rp and (query) point q, the maximum 
possible inner product between the point q and the ball Bp^ is bounded from above by: 

max {q,p) < {q,po) + Rp \\q\\ . (4) 

Proof. Suppose that p* is the best possible match in the ball Bp^ for the query q and rp be the Euclidean 
distance between the ball center po and p* (by definition, rp < Rp). Let 9p be the angle between the vector 
Po and the vector pop* , (j) and ojp be the angles made at the origin between the vector po and vectors q and 
p* respectively (see figure [SJ. The length of p* in terms of po and 9p is: 



b*ll = ^J{\\po\\+rpcos0py + {rpSmOpy. (5) 
The angle ujp can be expressed in terms of po and 9p as: 

IIpoII +''pC0s6lp . rpainOp 
= --n., = -j-^. (6) 

Let Oq^p* be the angle between the vectors q and p* . With the triangle inequality of angles, we have: 

\Oq,p* \ > 10- Wp|. 

Assuming that the angles lie in the range [— tt, tt] (instead of the usual [0,27r]), and the fact that cos(0) — 
cos{~9), we get: 

cos9q^p» < cos{(j) — ojp), (7) 
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since cos(-) is monotonically decreasing in the range [0,7r]. Using this inequahty we obtain the following 
bound for the highest possible inner-product between q and any point in the ball: 

max {q,p) = (<Z,P*)(by assumption) 

= ||g|| ||p*|| cos^^^p. 

< ||g|| ||p*|| cos((/- -Wp), 

where the last inequality follows from equation [T] Substituting equations [S] & [B] in the above inequality, we 
have 

max < \\q\\ {cos <j){\\po\\ + rp cos 9p) + sin (l){rpSm9p)) 



< 



\\q\\ max (cos '/>(||po|| + fp cos Op) + sin (/){rp sin 6'p)) 



= llgll (cos 0(||po|| + fp cos (p) + sin(j)(rp sin 4>)) 

< (cos0(||po|| + ^pCOS(/)) + sin0(i?psin(/))) (since rp < _Rp) 



The second inequality comes from the definition of maximum. The following equality comes from maximizing 
over 9p. This gives us the optimal value of 9p ~ (p. Simplifying the final inequality gives us equation S] □ 

For the tree-search algorithm (Alg. |4]), we set the maximum possible inner-product between q and a tree 
node T as 

MlP{q,T) = {q,T.^i) +T.R\\q\\ . 

This upper bound can be computed in almost the same time required for a single inner-product (since the 
norms of the queries can be pre-computed before searching the tree). This algorithm is evaluated against 
the naive linear search algorithm in section [S] 



4 Dual-tree based Search 

For a set of queries, the tree can be traversed separately for each query. However, if the set of queries is 
very large, a common technique to improve efficiency of querying is to index the queries in the form of a 
tree as well. The search is then subsequently done by traversing both trees simultaneously using the dual- 
tree algorithm [15 . The basic idea is to amortize the cost of tree-traversal for a set of queries which are 
very similar to each other and would follow (approximately) the same path down the tree. The dual-tree 
algorithms have been applied to different tree-based algorithms like nearest-neighbor search '15' and kernel 
density estimation [16 with provable theoretical runtime bounds |31) . 

4.1 Dual-tree Branch-and-bound Algorithm 

The generic dual-tree algorithm is presented in Algorithm [51 Similar to the Algorithm 31 the algorithm 
traverses down the tree on the reference set S (referred to as the RTree). However, the algorithm also 
traverses down the tree on the set V of queries (QTree), resulting in a four- way recursion. At each step, 
the algorithm is at a QTree node Q and a RTree node T. For every Q, the value Q.X denotes the minimum 
inner-product between any query in Q and its current best-match candidate. If this value is greater than 
the maximum possible inner product, MIP((5,r), between any query in Q and any reference point in T, 
this part of the recursion is no longer explored. When the algorithm is at the leaf level of both the trees, it 
obtains the best-matches for each query in the QTree leaf by doing a linear scan over the RTree leaf. 

In this section, we explore two ways of indexing the queries - (1) indexing the queries using the ball-tree 
(2) indexing the queries using a novel data structure, the cone-tree. In the following subsections, we derive 
expressions for MIP(Q,T) each of these kinds of QTree. 
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Algorithm 6 DualSearch(QTree Node Q, RTree Node T) 
if Q.A < MIP((3,r) then 

/ / This node has potential 

if tsLeafiT) & tsLeaf{Q) then 

for each q G Q.S do 
LinearSearch(g, T.S) 

end for 

Q.A <— miiiggQ.s g.A 
else if isLeaf(T) then 

DualSearch(Q.lc, T); DualSearch(g.rc, T); 

Q.A min{Q.lc.A, Q.rc.A} 
else if isLeaf{Q) then 

Ii ^ MIP(Q,r.lc); Ir ^ MIP(Q,r.rc); 

if /; < Ir then 

DualSearch(Q, T.rc); DualSearch(Q, T.lc); 

else 

DualSearch(Q, T.lc); DualSearch(Q, T.rc); 
end if 
else 

/ / best depth first traversal 

Ii 4- MIP(Q.lc,r.lc); Ir ^ MIP(Q.lc,r.rc); 

if Ii < Ir then 

DualSearch(Q.lc, T.rc); DualSearch(Q.lc, T.lc); 
else 

DualSearch(Q.lc, T.lc); DualSearch(Q.lc, T.rc); 
end if 

Ii ^ MIP(Q.rc,r.lc); /,. ^ MIP(Q.rc, T.rc); 
if Ii < Ir then 

DualSearch(Q.rc, T.rc); DualSearch(Q.rc, T.lc); 
else 

DualSearch(Q.rc, T.lc); DualSearch(Q.rc, T.rc); 
end if 

Q.A min{Q.lc.A, Q.rc.A} 
end if 
end if 

// Else the node is pruned from computation 



Algorithm 7 FindExactMaxIPDualTree(Query Set V, Reference Set S) 

T ^ MakeBallTree(S') 
Q ■(- MakeQueryTree(V') 

V trees nodes Q' in the tree Q, Q'.A < oo; 

V queries q £ V, q.hm 0, q.X i oo; 

DualSearch(Q, T); 

V queries q £ V , return q.bm; 



Figure 6: Dual-tree Search; The tree-building subroutine for the set of queries MakeQueryTree" can 
be the "MakeBallTree" subroutine (Alg. [2|) or the " MakeConeTree" subroutine (Alg. [9|). The object q.hm 
contains the current best-match for the query q. Q.X denotes the lowest affinity between any query in the 
node Q and its current best-match. The function MIP((5,T) denotes the upper bound on the maximum 
possible inner-product between any query in the node Q and any point in the node T. 



10 



X 



O 



Figure 7: Bounding between two balls. 
4.2 Ball Tree for Queries 

Theorem 4.1. Given two balls Bpg and Bq^^ centered at po o-nd with radius Rp and Rq respectively, the 
maximum possible inner-product with any pair of points p e Bp^ and q € S^" is bounded from above by: 

max {q,p) < {qo,po) + RqRp + \\qo\\ Rp + \\Po\\ Rq- (8) 

Proof. Consider the pair of point {p*,q*),p* S Bp^ ,q* e -B^" be such that 

= max {q,p). (9) 

Let 6p be the angle po makes with the vector pop* , and 9q be the corresponding angle in the query ball. 
Let Wp be the angle between the vectors po and p* and uiq be the angle between the vectors go and q* . Let 
rp be the distance between po and p* , Vq be the distance between qo and q* . Finally, let (p be the angle made 
between po and qo at the origin. 

Some facts for the ball Bp^ (the facts are analogous for tlw^ ball S,^'): 



Ibll = \/\\Pof + rl + 2\\po\\rpCosep, 



Poll +rp cos 9p . rpSinOp 
sma;„ — 



' Ibll ' ' \\P* 

Using the triangle inequality of the angles, we know that: 

\0q',p' \ > 10- {0Jp+UJq)\, 

giving us the following: 



{q*,P*) = Ibll Ik* II cos(<^ - {Up+ujq)) (10) 
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Replacing ujp and Uq with 9p and 6q by using the aforementioned equalities (similar to the techniques in 
proof for theorem 13. ip . we have: 



{q*,P*) = ((7o,Po) +rpr, cos((/)- (^p + 61,)) +rp llgoll cos( 



Tq llpoll cos( 



< 



{<10,Pa) + rpTq COs{(f) - {Op + 9q)) + Tp llqoll COS(0 - Op) + Tq \\po\\ 



< ma,x{qQ,po) + rpVq + rq\\po\\ + rp\\qo\\ (since cos(-) < 1), 

< {qo,PQ) + RpRq + Rq \\pa\\ + Rp \\qa\\ , 



(11) 

(12) 



where the first inequality comes from the definition of max and the final inequality comes from the fact that 



rp < Rp, Tq < R, 



□ 



For the dual-tree search algorithm (Alg. 15]), the maximum-possible inner-product between two tree nodes 
Q and T is set as 

MIP(g, T) = (qo,Po) + RpRq + Rq IIpoII + Rp \\qo\\ . 



It is interesting to note that this upper bound bound reduces to the bound in theorem 13.11 when the ball 
containing the queries is reduced to a single point, implying Rq = 0. 




Figure 8: Cone-tree: These cones are open cones and only the angle made at the origin with the axis of 
the cone is bounded for every point in the cone. The norms of the queries are not bounded at all. 



4.3 Cone-trees for Queries 

An interesting fact is that in equation [TJ the point p, where the maximum is achieved, is independent of the 
norm ||q|| of the query q. Let 9q r be the angle between the q and r at the origin, then the task of searching 
for the maximum inner-product is equivalent to search for a point p € S such that: 

p = argmax ||r|| cosflg^r- (13) 

This implies that we only care about the direction of the queries irrespective of their norms. For this reason, 
we propose the indexing of the queries on the basis of their direction (from the origin) to form a cone- 
tree (figure [8]). The queries are hierarchically indexed as (possibly overlapping) open cones. Each cone is 
represented by a vector, which corresponds to its axis, and an angle, which corresponds to the maximum 
angle made by any point within the cone with the axis at the origin. 
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Algorithm 8 MakeConeTreeSplit(Data Q) 
Pick a random point x € Q 
A -f- arg min^'es cos 0,^ ^/ 
B arg minx' gs cos 6'^_x' 
return {A,B). 



Algorithm 9 MakeConeTree(Set of items S) 
Input - Set S 
Output - Tree T 
T.S ^ S 
T.^ <— mean(S') 
T.C <~ miupes cos St.m.p 
if IS"! < A^o then 

return T 
else 

iA,B) ^ MakeConeTreeSplit(S) 
Si {p e S : COS 9 A,p > cos6b.p} 

Sr^S\Sl 

T.lc ^ MakeConeTree(5'i) 
T.rc MakeConeTree(5'r-) 
return T 
end if 



Figure 9: Cone-tree Construction: The object T.S denotes the set of points in the node T, T./x denotes 
the Euclidean mean of the items in the node T and T.C denotes the cosine of the maximum angle made by 
any point in the node with T./x at the origin. The angle made between any two points A and B at the origin 
is denoted by 9a,b- 

4.3.1 Cone-tree Construction 

The cone-tree construction is very similar to the ball-tree construction. The only difference is the use of 
cosine similarity instead of the Euclidean distances for the task of splitting. The cone-tree construction 
pseudo-code is presented in Figure [H 

4.3.2 Cone-Ball Bound 

Since the norms of the queries do not affect the solution in equation I13[ we assume that the norms of the 
queries are all equal to 1 for convenience. 

Theorem 4.2. Given a ball Bpg of points centered at po with radius Rp and a cone Cq^ with the axis of the 
cone go o^i^d apertur^ of 2LUq > 0, the maximum possible inner-product between any pair of points p € Bp^ , 
q S Cqg is bounded from above by: 

max {q,p) = max ||p||cos6'qp 

< \\po\\cos{{\c^\-ujq}+) + Rp, (14) 
where 4> is the angle made between po and qo at the origin and the function {x}+ — maxjo;, 0}. 
Proof. There are two cases to consider here: 

(i) \(l)\<UJq 

^The aperture of the cone is twice the angle made between the axis and the perimeter of the cone. 
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y 




Figure 10: Bounding between a ball and a cone 



(ii) \4>\>ujg 

For case (i), the center po of the ball Bp^ lies within the cone Cq^ ^ implying that 

max \\p\\cofiOq,p <\\po\\ + Rp. (15) 

since there could be some query q* £ Cqg which is in the same direction as po, giving the maximum possible 
inner-product. 

For case (ii), let us assume that > without loss of generality. Then (f> > ujq. Continuing with the 
similar notation as in theorem 13 . 1 1 fc HTT] for the best pair of points {q* ,p*) as well as the notation from figure 
[TOl we can say that 

\0p'.q- \ > If/) - - UJp\ (16) 

Since coq is fixed, we can say that 

max ||p||cos6'gp < ||p*|| cos6'g. p- (by def.) 

< \\p*\\cOs{(f)- UJq- UJp). (17) 

Expressing ||p*|| and LUp in terms of ||po|| jfp and 6p, and then subsequently maximizing over 6p and using 
the fact that rp < Rp, we get that 

max IIpII cos 6',,p < IIpoII cos(0 - Wg) + i?p. (18) 
Combining case (i) and (ii), we obtain equation 1141 □ 



5 Experiments and Results 

In this section, we evaluate the efficiency of the two proposed algorithms [S] & [71 For the dual-tree algorithm, 
we use the two variations - (i) the set of queries indexed as a ball-tree (referred to as Alg. El^B)), (ii) the 



14 



Dataset 


© 




1^1 


Bio 


74 


210,409 


75,000 


Corel 


32 


27,749 


10,000 


Covertype 


55 


431,012 


150,000 


LCDM 


3 


10,777,216 


6,000,000 


LiveJournal 


25,327 


121,625 


100,000 


MNIST 


786 


60,000 


10,000 


MovieLens 


51 


3,706 


6,040 


Netflix 


51 


17,770 


480,189 


OptDigits 


64 


1,347 


450 


Pall7 


7 


100,841 


100,841 


Physics 


78 


112,500 


37,500 


PSF 


2 


3,056,092 


3,056,092 


SJ2 


2 


50,000 


50,000 


U-Random 


20 


700,000 


300,000 


Y!-Music 


51 


624,961 


1,000,990 



Table 1: Datasets used for evaluation: The dimensionality 'D and number of points in the reference set 
S and the set of queries V. 

set of queries indexed as a cone-tree (referred to as Alg. EJC)). Since we are not aware of any efficient 
exact method for maximum inner-product search, we compare our proposed algorithms to the linear search 
algorithm (Alg. [S]). We report the speedup of the proposed algorithms over hnear search. Speedup is defined 
as the ratio of the time taken by the linear search and the time taken by the evaluated algorithm. For the 
trees, the leaf size Nq can be selected by cross-validation (choosing the leaf size giving the highest speedup). 
However, for our experiments, we choose a ad hoc value of Nq — 20 for all datasets to demonstrate the gain 
in efficiency without any expensive cross-validation. 

Datasets. We use a variety of datasets from different fields of data mining. We use the following collab- 
orative filtering datasets: MovieLens [17], Netflix [3] and the Yahoo! Music [12 datasets. After the matrix 
factorization stage, the matrix of item- vectors is used as the reference set and the matrix of user-vectors is 
used as the set of queries. For text data, we use the LiveJournal blog moods data set (19, . We also use 
the MNIST digits dataset [21] for evaluation. We also use three astronomy datasets - LCDM [27], PSF 
and SJ2. A synthetic data set (U-Rand) of uniformly random points in 20 dimensions is used to evaluate 
the performance of the tree-based algorithms on data sets without any underlying structure. The rest of 
the datasets are widely used machine learning data sets from the UCI machine learning repository |; 5 . The 
details of the datasets are presented in Table [Tj and the size of the datasets (in bytes) is presented in figure 
[TTj For the collaborative filtering datasets, there is a clear definition of the reference set (the items) and 
the set of queries (the users). For the rest of the data sets, we randomly split the datasets into query and 
reference sets. 

Tree Construction Times. The tree-building procedure is extremely efficient. We present the tree 
construction times in table [5] and contrast them with the runtime of the linear search algorithm (Alg. [3]). 
For some of the larger data sets, the extrapolated runtime of Alg. |3| is reported. In the last column, we 
present the ratio of the tree construction times with the runtimes of Alg. |3l For algorithm [5] & |71[B), the 
tree construction involves building one and two ball-trees respectively. For algorithm |7]JC), the queries are 
normalized to have unit length for convenience since the norms of the queries do not affect the answers 
(equation I13p. Following the query normalization, two trees are built. We include the query normalization 
in the tree construction time for completeness. This is the reason for the significant difference between 
construction times for algorithm |71[B) and|7l[C). 

The numbers in the last column of table |5|(R) show how small the construction times are with respect to 
the actual linear search. The highest ratio is 0.15 for the OptDigits dataset. This implies that any speedup 
over 1.18 at search time is enough to compensate for the tree construction time. For most of the datasets, 
this ratio is much lower. Moreover, this tree building cost is a one time cost. Once the tree is built, it can 
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Figure 11: Total dataset sizes (in bytes): The combined sizes of the reference set and the query set for 
each data set are presented in this figure. 

be used for searching the dataset muhiple times. 

Search Efficiency. The speedups over hnear search are presented in Table \5\ We have reported every 
dataset we evaluated our algorithm on. Overall, the speedup numbers vary from as low as 1.13 for the 
OptDigits dataset to over 10^ (4 orders of magnitude) for the LCDM and the PSF dataset. An important 
thing to note here is that for datasets with low speedup (below an order of magnitude) with Alg. \E\ the 
speedup numbers for all three algorithms were pretty low and fairly comparable for all three algorithms. 
However, even a speedup of 2 is pretty significant in terms of absolute times. For example, for the Yahoo! 
music dataset, a search speedup of mere 2 with a tree construction time of 120 seconds gives a saving of 19 
hours of computation time. For most datasets with a high value of speedup for Alg. [Sj the speedups for the 
dual-tree algorithms are also very high. 

There are three important things to note here. Firstly, the dual-tree algorithms (Alg. [7]) do not perform 
very well if the single-tree algorithms (Alg. [S]) does not have a high speedup. This is mostly because the tree 
is unable to find tight bounds and hence has to travel every branch. The dual-tree scheme loosens the bound 
to amortize the traversal cost over multiple queries. But if the bounds are bad for algorithm [Sj the bounds 
for the dual-tree are much worse. Hence, the dual-tree algorithm does not show any significant speedup. 
Secondly, the dual-tree algorithm (especially Alg. EJC)) starts outperforming the single-tree algorithm 
significantly when the set of queries is really large. This is a usual behavior for dual-tree algorithms. The 
query set has to be large enough for the gains from the amortization of query traversal of the reference tree 
(RTree) to outweigh the computational cost of traversing the query-tree (QTree) itself. Finally, the dual-tree 
algorithm with ball-trees for the query set is generally significantly slower than the dual-tree with a cone-tree 
for the queries. There are possibly two possible reasons for that - (i) The cones provide a tighter indexing 
of the queries than balls. A single cone can be used to index points in multiple balls which lie in the same 
direction but have varying norms, (ii) The upper bound for MIP(Q, T) in equation [TT] is fairly loose. We do 
provide two ways of obtaining tighter bounds in the Appendix, but we have not yet evaluated the algorithm 
with the new bounding techniques. 

We also consider the general problem of obtaining the points in the set S with the k highest inner-product 
with the query q. This is analogous to the fc-nearest neighbor search problem. We present the speedups of 
our algorithms over linear search for fc = 1, 2, 5 & 10 in figure [T^ 
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Dataset 


Algig 


AlglZlB) 


AlglUC) 


Alg|3] 


R(%) 


Bio 


4.3 


5.7 


10.6 


4,028 


0.25 


Corel 


0.2 


0.27 


0.66 


43 


1.5 


Covertype 


5.5 


7.2 


14.8 


14,885 


0.1 


LCDM 


36.7 


56.46 


99.3 


1,984,200 


0.005 


LiveJournal 


2223 


4073 


4745 


517,194 


0.92 


MNIST 


8.06 


9.1 


11.38 


817 


1.5 


MovieLens 


0.03 


0.08 


0.27 


4.62 


6 


Netflix 


0.2 


8.27 


33.5 


1,878 


1.7 


OptDigits 


0.01 


0.012 


0.022 


0.135 


15 


Pall7 


0.26 


0.52 


1.4 


364 


0.4 


Physics 


2.33 


3.0 


5.8 


1,114 


0.5 


PSF 


9.06 


18.1 


34.95 


282,514 


0.01 


SJ2 


0.1 


0.2 


0.46 


75 


0.6 


U-Rand 


4.94 


6.9 


15.64 


26,586 


0.6 


Y! Music 


9.72 


28.85 


112.5 


137,306 


0.08 



Table 2: Tree construction time (in seconds) contrasted with the linear search time (in seconds). 



Dataset 


AlgE] 


AlglUB) 


AlgiaC) 


Bio 


7,059.62 


6.55 


273.52 


Corel 


14.27 


17.38 


7.68 


Covertype 


927.51 


10.05 


773.34 


LCDM 


29,526 


1,327 


101,950 


LiveJournal 


28.04 


10.42 


15.45 


MNIST 


2.61 


2.22 


2.5 


MovieLens 


2.23 


1.36 


1.67 


Netflix 


1.98 


1.92 


1.84 


OptDigits 


1.13 


1.10 


1.10 


Pall7 


1,020 


23.14 


2,285 


Physics 


4.93 


4.0 


4.08 


PSF 


61,502 


96,570 


125,800 


SJ2 


544 


190 


767 


U-Rand 


3.76 


3.18 


3.28 


Y!-Music 


2.11 


2.09 


2.16 



Table 3: Speedups over linear search for k ~ 1. 



6 Max-kernel Operation with General Kernel Functions 

In this section, we show a method to apply the proposed algorithms in a inner-produet space where the 
inner-products are defined by a kernel function, but it is not possible to explicitly represent the points in 
the (^-space. 

Without an explicit representation, the tree construction has to be modified since there would be no 
explicit representation of the mean of a set. For a tree node T with the set of point T.S, the mean in 95-space 
is defined as 

1 



pGT.S 

II might not have an explicit representation, but it is possible to compute inner products with /i as follows: 

However, this computation is possibly very expensive (as opposed to the operation in equation |4] which is 
equivalent to a single inner-product). Instead of picking the mean of the set in the ip-space as the center of 
the ball, we propose picking the point in the <p-space which is closest to the mean /x as the new center. So 
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Figure 12: Speedups over linear search for k = 1, 2, 5 & 10. 



>- 



the new center pc is given by: 



Pc = arg^min^||(^(r) - ^1 



2 



2 

^''^r^TS^'-'''^^ ~ IT^ H /C(r',r). (19) 

' ■ ' r'eT.S 



This operation is quadratic in computation time, but is done at the preprocessing phase to provide efficiency 
during the search phase. Given this new center pc, we can compute the radius Rp of the ball enclosing the 
set T.S as follows: 

RI = max ||v3(r) - (^(pc)ll^ 

= maxJC{pc,Pc)+IC{r,r)~2IC{r,pc). (20) 

r^T.S 

Now given this method of choosing the center and evaluating the radius, a ball-tree can be built in 
the (/3-space using Algorithm [2] without ever requiring the explicit representation of the points. Given this 
ball- tree, the equation |4] in theorem 13. II can be modified to this situation as follows: 



MIPiq,T)=ICiq,p,)+Rp^lC{g;qj, (21) 
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where Pc is defined in equation [TO] and Rp is defined in equation [501 Computing this upper bound is equivalent 
to a single kernel function evaluation {lC{q, q) be pre-computed before searching the tree). Using this upper 
bound, the tree-search algorithm (Alg. [S|) can be performed in (^-space without any explicit representation 
of the points. We will present the evaluation of this method in the longer version of the paper. 

Using the same principles, the dual-tree algorithm (Alg. [7]) can also be applied to the tp-space without 
any explicit representation of the points. For the dual-tree with ball-tree for the queries, the upper bound 
on the maximum inner-product between queries in node Q and points in node T in theorem 14. II becomes: 



MIP(Q,r) = K.{q,,p,) + RpRq + Rp^lC{q,,q,) + Rq^}C{p^,p,), (22) 

where Pc and qc are the chosen ball centers in the 95-space with radius Rp and Rq respectively. 

For queries indexed in a cone-tree, the central axis of the cone can be the point in the (^-space making 
the smallest angle with the mean of the set in the iy9-space. Since the queries are supposed to be normalized 
in the (^-space, for a query tree node Q, the mean of the set Q.S is supposed to be: 



So the new central axis qc of the cone is given by: 



arg max 



q€Q.S \\fl\\ \\ip{q)\\ 



q'&Q.S ^'^il'^l') 

= arg max , . (23) 

Again, this computation is quadratic in the size of the dataset, but provides efficiency during search time. 
The cosine of half the aperture of the cone is now given by: 

{'p{Qc),(p{q)) 

cos ujn = mm 



qeQ.S \\ip{qc)\\ \\fiq)\\ 



leQ.S ^IC{qc,qc)ICiq,q) 



(24) 

Given qc and ujq, the upper bound in theorem 14.21 for a cone-tree node Q of queries and a ball-tree node T 
of reference points is given by: 

MIP(g, T) = ^ICipc,Pc)cos{M - ujq}+) + Rp, (25) 

where (p is defined as: 

JC{pc,qc) 

cosr — 



y^JC{qc,qc)IC{pc,Pc) 

This bound is very efficient to compute as it only requires a single kernel function evaluation (the terms 
JC(j>c,Pc) and JC{qc, qc) can be pre-computed and stored in the trees). 

7 Conclusion 

We consider the general problem of maximum inner-product search and present three novel methods to 
solve this problem efficiently. We use the tree data structure and present a branch-and-bound algorithm for 
maximum inner-product search. We further extend it to the case where the set of queries is very large. We 
evaluate the proposed algorithms with a variety of datasets and exhibit their computational efficiency. 

A theoretical analyses of these proposed algorithms would give us a better understanding of the compu- 
tational efficiency of these algorithms. We do not have any rigorous runtime bounds for our algorithm and 
it would be part of our future work. 
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A Tighter Bounds with Optimization 

In this section, we present two ways to get a tighter bound on equation [TT] with respect to 9p and 9q. The 
maximum inner product bound MTP{Q,T) between two balls is given as: 

{q*,P*)< max {qo,po) +rprqCOs{(f>~ {9p + 9q)) + rp\\qo\\cos{(l)- 9p) + rq\\po\\cos{(f>~ 9q). (26) 

A.l Two- variable Optimization 

Assuming that 



{q*,P*) < ina.x{qo,po) + RpRqCOs{(j) - {9p + 9q)) + Rp\\qo\\cos{<j) ~ Op) + Rq\\po\\cos{4> - 9q) (27) 



\cl^-{0p + 9q)\<-,\<t>~9p\<-,\^^9q\<-, 



we can say that: 




(28) 



Now^^^ 




sin(0 - Op) Rq 



(29) 



sin(</. - {9p + 9q)) 
sm{(j) - 9q) 



Iboll 



(30) 
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e,) 


del 




d'fiep, 


e,) 


del 




d^ie^, 


e,) 



The second order conditions are the following: 

= -RpRqcos{(f)-{ep + eq))- Rp\\qo\\co&{(t>^ep), (31) 

= -RpRqCOs{(j)-{ep^eq))-Rq\\pQ\\cOs{(j)-eq), (32) 

^^^^^^ = -RpR,cos{cj)-{ep + e,)), (33) 

which are all < for the stated range of ^, Op and 9q . So the optimal values obtained from the optimality 
conditions (equations [29] & 1301) correspond to the maximum. However, the optimality conditions do not 
have an analytic solution for Op and Oq. Hence, any efficient optimization algorithm can be used to solve 
maxgp^g^j f{ep,9q) in the specified range. 

A. 2 One-variable Optimization 

Another approach is the following: 

max (po,(7o> + rprqCos{cj) - {dp + Oq)) + rp| |go| | cos((/) - Op) + rg||po|| cos((/) - Oq) 

T (34) 
< max maxpp go + rpVq cos{(p - {Op + Oq)) + rp||go|| cos(^ - Op) + r,||po|| cos(^ - Oq), 



^q,' p,' q ^p 



Since for fixed Oq, ujq is fixed. And for a fixed Uq, using the single-tree bounding. Op — {(p — ujq). Making this 
substitution in equation 1341 we get the following optimization task: 

{p',q*) < T^Si-x {po,qo) +rp\\q*\\+rq\\po\\cos{(l)-Oq) (35) 

< niax(po,go) + ^glboll cos((/) - Oq) + Rp^JWqoW^ + + '^Rq\\qo\\ cosOq (36) 
= fiOp), 

where the second inequality comes from the assumption that 

1^,1 <f, 10-0,1 <|, 

and the fact that rp < Rp, Vq < Rq. 

The first-order optimality condition gives us the following: 

RqWqoWsmOp II • / , a \ 

Rp I = Rq\\po\\?'\n{(j)- Oq), 

^Il9o|p+i?2 + 2i?,||(zo||cos^, 
while the second-order derivative is given by: 

DON i|Cos6'p(||go|P+i?g + ^<;lko||cos6',)-Ki?,||go|l „ 11 ,, „, 
-RqRpWqoW ^ -.ry^ h i?, | |po 1 1 cos(0 - 6*^, 

{\\qo\\^ + Rl + 2Rq\\q^\\cos0qf^ 

which is always < implying that the optimal Op is the maximum even though the equation IA.2I does not 
give an analytic solution for Op. An efficient optimization algorithm can be used to solve this one dimensional 
optimization problem maxg^ f{Op) to obtain tight bounds for MIP(Q, T). 
promising result. 
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